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ABSTRACT 
It is relatively easy to generate, by digital computer, 
large numbers of seemingly independent random numbers with a 
uniform distribution over a fixed range, say - y < X < y. 
Methods of generating gaussian, or normal, random numbers generally 
are based on either non-linear transformations on random numbers 
from a uniform population, or the summing of enough independent 
numbers from a uniform population for the central limit theorem 
to be applicable. In the first case a time-consuming evaluation 
of a complicated function is involved. The second method is also 
slow because a large number of uniform random variables must be 
generated and summed for each normal random variable obtained. 
This note discloses a method based on the central limit theorem, 
except that the summing of N uniform random variables gives N 
normal random variables. The approach is to form an N dimensional 
vector whose components are uniform random variables, multiply the 
vector by a Hadamard matrix, and use the resulting components as 
normal random variables. It can be shown that the resulting N 
components have a uniform density inside a N-dimensional hypercube 
aligned with diagonals along the coordinate axes. However, the 
one dimensional marginal densities, the two dimensional marginal 
densities, indeed all the marginal densities tend toward the normal 
density as N gets large. Furthermore the components are uncorre- 
cted and have equal variance, independent of N. However, some 
of the fourth moments, which should be zero for independent 
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normal random variables, are not zero for our derived set 
(although these moments do approach zero as N becomes large) . 
These moments can be made zero, however, by randomly changing, 
or not changing, the sign of each component. 

The method proposed is very fast because the principal 
step, Hadamard matrix multiplication, requires only N log2 N 
additions to produce N components. 



Accepted for the Air Force 

Franklin C. Hudson 

Chief, Lincoln Laboratory Office 



IV 



A NEW METHOD OF GENERATING GAUSSIAN RANDOM 
VARIABLES BY COMPUTER 

Many scientific computations require the generation by a 
computer of a large number of seemingly random numbers with a 
multidimensional normal (gaussian) probability density function. 
However, although there are efficient algorithms available for 
generating uniformly distributed random numbers, the algorithms 
which generate normal random numbers are generally much slower. 
One typical algorithm begins by generating a moderate number of 
uniform random variables (say N) and forming the sum, which, for 
large N, is approximately normal, according to the central limit 
theorem. This method is slow because many uniformly distributed 
random variables must be developed for each resulting normal 
random variable, and because many additions must be performed 
for each normal random variable, as well. Another typical method, 
more attractive because it is less dependent on approximation, 
begins with a single uniform random variable and performs a non- 
linear transformation of the random variable to form a new random 
variable with a normal density. Unfortunately, the non-linear 
function which must be evaluated is complicated and cannot be 
computed quickly on most digital computers. 

This note describes a variation of the central limit theorem 
approach, in which N uniform random variables (r.v.) are trans- 
formed into N approximately normal r.v. by a Hadamard transformation, 

T~. The term "pseudorandom" is often used to describe the determin- 
istically generated, but seemingly random numbers. 



This operation takes y N log 2 N additions and y N log 2 N subtractions 
which is to say, log ? N additions or subtractions per normal r.v. 
obtained, and, of course, it is only necessary to generate one 
uniform r.v. for each normal r.v. obtained. The plan of the 
paper is to first describe the ideal uniform and normal multivariate 
densities, second to describe the Hadamard transformation used, 
third to derive the properties of the transformed variables, and 
fourth, to consider the properties of the transformed variables 
after they have been subjected to random sign changes. 

I. Properties of the Ideal Uniform and Normal Densities 

We define the ideal uniform r.v., X , as one whose probability 
density is 

1 1 l*„l < I 

( ° l X „l » 7 
For several such r.v.s X Q ,Xj, . . . »X N i to be independent implies 
that the joint density is 

I 1 if all |X n | < \ 

P xx...x (X 0' X l 5X 2 5 *** 5X N-l ) ] 1 (2) 

if any |xj > j 

On an N dimensional space, the joint density is unity on the 
inside and zero on the outside of a hypercube centered on the 
origin with coordinate axes passing from the center of each face 
to the center of the opposite face. 



The expected value, E(X ) = X , is zero and, by integration, 

all the other moments can be found. For the case of even order 

moments of X , 
n' 

?P = -*— i (3) 

n 2 zp (2p+l) 

and for odd order moments 



X 2p+1 = (4) 

n v ' 



Since the mean is zero, the concept of central moments is super- 
fluous. For joint moments, we appeal to independence, which says 
that the moment of a product is the product of the moments. 
Therefore 



(xX X c> " 0£)(x£)Otc> (5) 

and, of course, these joint moments are zero except when m, n, I 
are all even. 

It is worth taking note of the second moment, X = 1/12. 

The normal probability density of an r.v. Y , with zero mean, 
is given by 

W - ~=^- e _ V 2a (6) 

y n /Zt? a 



2 
where a is the variance. We shall be comparing several densities 

in this note and we shall force the variances to all be equal to 

1/12. Therefore a 2 = 1/12. 

The multivariate density of several such independent variables, 

Y ,Y 1' ' ' ' ' Y N-1' 1S 



P yyy..y (Y 0» Y l''"' Y N-l ) 

1 -(Y2 + Y 2 1+ ... + Y2_ 1 )/2a 2 

TZTF e (7) 



The odd order one dimensional moments of Y are 

n 



7^= o 

n 



(8) 



while the even order moments may be shown to be 

^2p = l'3'5'"(2p-l) . l-3-5---(2p-l) (g) 

(l/a) 2p 12 p 

Again, by independence, the joint moments are equal to the product 
of the individual moments, so that the moment 



#J = o£)(y{J)(y*) (10) 



is zero unless m, n, I are all even. 



The normal density is never zero, and is, in fact, radially 
symmetrical, i.e., P (Y n ,Y-. , . . . jY^ i ) is only a function of 
the distance of Yq,Y-j, . . . ,Y N i from the origin in N-space. This 
is also a true statement for any marginal density. The same 
statement is not true for the multivariate uniform density, which 
takes on different probabilities in the directions of the vertices 
of the hypercube than along the direction of the faces of the cube. 

II. The Hadamard Matrix 

By definition, a Hadamard matrix is an orthogonal matrix 
whose elements are all either +1 or -1. Hadamard matrices do not 
exist for all dimensionalities; for example there is no 3 x 3 
Hadamard matrix. It is relatively easy to construct an N X N 
Hadamard matrix, however, if N is a power of two. We shall content 
ourselves with a particular Hadamard matrix for each power of 
two, for our purpose is not to investigate Hadamard matrices but 
to apply them. Letting N = 2 , and numbering our dimensions in 
the binary system 



i = i Q + 2i x + 4i 2 + ... + 2 c " 1 i c _ 1 * i * N-l 



j = j Q + 2j x + 4j 2 + ... + 2 c - 1 j c _ 1 * j * N-l 



the elements of our particular Hadamard matrix, H, are given by 



(id 



h . (.d^o ^hh ^hh ... ^Wc-i 



(12) 



We shall show in an appendix that multiplication of an N component 
vector by H can be accomplished with N log 2 N additions and 
subtractions . 

It is clear that h. . = h. . so that H is a symmetric matrix. 
Also, the zero row and the zero column are composed of +l's 
only. The product of the elements in the i row with the k row 
give the elements of another row of the Hadamard matrix. 



h. . h, , - h f . (13) 



where i is the row whose number is 



l = (i Q © k Q ) + 2(i x 9 k x ) + ... 



+ 2 C " 1 (i c _ 1 k^) (14) 



where 9 means "exclusive or". We shall say as a shorthand, 

I = i © k 

in place of (14). By symmetry, the product of columns is also a 
column. 



A consequence of (13) is the orthogonality of the Hadamard 
matrix. Strictly speaking the matrix 



7n h 



is not only its own transpose but its own inverse, and it is the 
matrix which we shall use in the transformation of uniform r.v. 
to normal r.v. 

We give as an example the 4x4 Hadamard matrix 



H 
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1-1 1-1 
1 1-1-1 
1-1-1 1 



III. Hadamard Generation of Almost Normal r.v. 

We begin with a set of uniformly distributed r.v. assumed 
independent. We form a new set of r.v. from the uniform set by 
the operation 



W m = 7N 



N-l 

S 
n=0 



X 



n,m n 



m = 0,1, . . .,N-1 



(15) 



and we shall be interested in the properties of the set of r.v. 
W m . For any single random variable, it is easy to see that P (W m ) 



m 



is the density of a sum of independent r.v. with zero mean and 



equal variance. It is known from the central limit theorem that 
such a sum approaches a gaussian r.v. in the limit of large N. 
The approximation is good in the region of W near zero, bad in 
the tails. It is, of course, clear that 



|W m l * \ M 



(16) 



so that P (W ) is zero outside this range. However, for large N, 
t v/N is many standard deviations from the mean, so that the 
difference can be tolerated. Indeed, the approximation is the 
same as for the common method of generating normal r.v. except 
that we hope to use a larger value of N and so obtain a good 
approximation. 

We have generated, however, N different r.v. from the same 
N uniform r.v. and we should not expect the W f s to be independent. 
This leads us to investigate the multivariate density 

^...wVi Vi)- But 

P w...w< w 0' W l'-"' W N-l> - P xx...x< X 0> X l""> X N-l>/ l°l 
where 



D 



aw. 



aw 



aX o ' ••" aX N-l 



a Vi aw N-i 



ax, 



ax. 



P ww...w< W 0' W l>"" W N-l> - P xx...x< X 0> X l>"-> X N-l> 



(17) 
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since the magnitude of the determinant of 777 H is unity. Of 
course, in (17), P (X n ,X-. , . . . ,X M J is to be evaluated at 



x 







X, 



L X 



N-l 



7ST 



H 



W 







w 1 



L w 



N-l' 



but this is only a change of coordinate system. In fact, there- 
fore, the N -dimensional joint density P is not normal but 
' J J ww. . ,w 

uniform. The density is one inside and zero outside a hypercube. 
This hypercube, however, is rotated. Whereas for the variables 
X , the hypercube was situated with the coordinate axes passing 
through its faces, for the W ' s we shall see that the coordinate 
axes pass through vertices of the hypercube. In fact, the 
coordinates of some of the vertices of the hypercube surrounding 

P xx...x (X 0» X l X N-1 ) are 



X. - ? 



i»J 



for any i, and such a vertex becomes rotated into the position 
with coordinates 



W. - 



/N i = i 



j ^ i 



which lies on the coordinate axis of the i coordinate in W-space. 
The vertex at -X- rotates to the opposite position, thus proving 
that the coordinate axes of W-space pass through diagonals of the 
hypercube. 

A hypercube in N -dimensions has 2 vertices. There are, 
however, only N coordinate axes, so only 2N vertices can lie on 
the coordinate axes. The remaining 2 - 2N vertices of the 
hypercube do not map onto the coordinate axes. At least some of 
the vertices are found in the "interior 11 of hyperquadrants , but 
not every hyperquadrant of the N-space can be lucky enough to have 
a vertex within it, for there are 2 hyperquadrants and at most 
2 - 2N vertices to go around. Therefore, we learn that 
P (Wq,W-,, . . . ,W N i ) is not only radially asymmetrical but is 
not symmetrical in each hyperquadrant of W-space. This shows up 
when we compute the moments of P . However, the most important 

WW . . .W ' 



moment of P w , . namely W W , is zero, showing that the variables 

ww . . »w 7 J m n' ' ° 

W are at least uncorrelated. We shall now compute all the first 

m r 

and second moments of P rTr , „(W n »W-, , . . . ,W M -.). 

WW...W X U 7 1 7 N-l' 

The first moments are all equal and are all zero, 



W n, " k S Vm X n " 75T Z V m X „ " ° < 18 > 



n ' n 
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because X =0. The second moments are 
n 



x N-l N-l 



WW, =i E Eh h K XX (19) 

a b N m =n „=n m > a n,b ^m n ' 



m=0 n=0 



. _2 



but XX is zero unless m = n, in which case it is a . Thus (19) 

m n ' ' 

becomes 



i N-l 9 

W W, = =r S h , mKN a (20) 

a b N m, (a©b) v ' 

m=0 

f-Vi 

By definition of a © b, it is zero if a = b, and the zero column 
of a Hadamard matrix consists of N ones. Therefore 



T 2 

W Z = a (21) 

a v ' 



as expected. But if a / b, a 6 b is not zero and sirce all the 
other columns of H contain as many -l's as l's, the sum 



£ h 
m 



m,(a©b) 



is zero, giving 



W a W b =0 (22) 



as stated above. 
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It is not difficult to argue that all the odd order moments 

of P TTr , TT are zero. We focus now on fourth moments. By definition 

WW . . .w J 

-i N-l N-l N-l N-l 
±7 Z Z Z Z 

N i=0 j=0 k=0 £=0 



W W,W W . = -tv z z z z 

abed ,,Z 



h i,a h jA,c h i,d¥jVl < 23 > 



To evaluate the fourth moments we note that X.X.X, X. = for 

i J k & 

most combinations of i,j,k,£ with the exceptions listed below; 



1. 


i=j, k«A, i#c 
i=k, j = A, i£j 
i=A, j=k, i^j 
i=j=k=A 


for which 
for which 
for which 
for which 


2 2 


a = 14T 


2. 


2 2 

x i x j - 


> - 1 
a " T4T 


3. 
4. 




^ 4 1 
a " T4T 

1 
" 80 



These are all disjoint cases. Therefore 

2/n-1 /N-l 



W a W b W c W d » 7 (T?) Lf h i,a9b k f Q h k,c9c 



x x 2/N-l /N-l 

+ ^ <l2"> lj h i,afcl jf Q h k,bO d , 
X Vi 

l l 2/n-1 /n-1 

+ 7 <I7> ^ h i,aedLf 9 h j,b$cj 

1 1 N_1 
+ ?W .f h i,a0b9c8d (24) 
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Equation (24) looks very complicated but it may be straightforwardly 

evaluated for any special case. Thus, for W we have a=b=c=d and 

a 

we find 



w|| - T" <TT) (n 2 " n ) + \ Cts) ( n2 - n ) + 



T 7T = 



ij- (i) (N 2 -N) + iy (i)N 

'I olio / ^ 

w a = wr + n ( m " wr) = 3a ■ ton (25) 



Similarly, if we let a=c , b=d, a^b we get 



2 



V? = 7 <ET> (n 2 -n) + iy (m)n 



XT 1 1 .4 a 2 



W aft - T4T + Tt0N - a + T5N < 26 > 

These results approach the gauss ian independent results 3a , a 
for large N. But consider the curious case of a,b,c,d all unequal 
but with aQbGcGd = 0. For this 



W a W b W c W d - SM < 27 > 
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An example of such a fourth moment is WqW-jW^Wo. This 
moment also approaches the expected result for independent normal 
r.v. as N - °°, However, it is somewhat disturbing that it is 
non-zero. We can, after all, convince ourselves that if a 
certain moment is within, say, 27o of the expected result for 
normal independent r.v. that this is somehow good enough. But 
what can we say if the expected moment should be zero? It is 
not clear how important this is, and indeed it probably depends 
on the application for which the r.v. are desired. 

Before considering the effect of a random sign change, it 
is well to summarize what we have shown so far. We have proposed 
a method by which uniformly distributed zero-mean random noise 
can be converted to approximately normal zero mean random noise. 
The one dimensional marginal densities are sums of N independent 
uniform r.v. and thus tend nicely toward normal for large N. 
The r.v. are also uncorrelated (in the sense that the expected 
value of a product is zero). Furthermore, based on the moments 
we have computed, the difference between a given joint moment 
and the same joint moment for independent gaussian noise seems 
to go to zero as 1/N. N can reasonably be made quite large since 
the amount of computation is proportional to logo N per normal 
r.v. generated. 
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IV. HADAMARD GENERATION OF MORE NEARLY INDEPENDENT RANDOM 
VARIABLES 

Let us now suppose that we have generated N random variables 



W . We have seen that some of the moments W W U W W , , which 
m abed' 

would be zero for an independent gauss ian set, are non-zero. We 
can form a set of N new random variables from the W ! s by changing 
the sign of each W, or not changing it, at random. That is, 
consider an easily obtained set of r.v., r , independent of each 
other and of W , with equal likelihood of being +1 or -1. For 
these r.v. the moments are all either zero or one. Even order 
moments are equal to 1, 



r r = 1 
m 



odd order moments are equal to zero, 



r 2 P + L - 

m 



and, by independence, the moment of the product of different r f s 
is the product of the moments. Now consider the set of random 
variables 



m mm 
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We shall compute some of the moments of the joint density 

P v^ V 0' V 1 , * ### ^N-1^ to s k° w that this density is nearly 
gaussian, for large N. First consider the moments which 



violated the gaussian assumption for the W f s. The V V, V V, is 
equal to 

r Q r K r ^ r ^i w Q W K W ^= (r Q r K r rj (W W K W W,) 
abcdabcd v abcd / v abcd / 



and unless a, b, c, d are equal at least in pairs, r r ,r r, will 
be zero. A similar argument will show that all joint moments of 

P which would be zero for independent, zero mean gaussian 

w • • • V 

noise are indeed zero. However, all other moments will be 



— 1 2 
unaffected by the conversion of W to V . Thus V = a , 

J m m a ' 



V V u = 0, and 
a b 



V aV - W aV - * 4 + T?N 



V a " - V * 3 * 4 " T§N 



We next consider the 6 order moments 



(29) 
2 



abcder abcder abcder 
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Discarding the trivially zero moments, we have for the other 



1 

N i j k i m n 



a b c d e f ^2 . . , i,a j,b Tc,c A,d 



h h £ X.X.X, X„X X 

m,e n,f ijkAmn 



We must again consider moments X.X.X^X-X X , most of which are 

xjKjfcmn' 

zero. The exceptions are listed below with the equal subscripts 
grouped and the ungrouped subscripts unequal. Thus 

Disjoint Cases Moment 

i j k l m n 2^ 



i 


j 


k 
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SL 
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i 


m 


k 
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m 
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.1 
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I 



1 
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Disjoint Cases Moment 
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This table allows us to compute all the 6 order moments 
for special cases. For example, if we were to compute V , we 
would have as part of the computation a sum 

1 \ l h 7^ - J^ i 

n 3 " i=o i ' (a • a ® a e a 9 a 9 a) l ' 448 i? 

and fifteen sums similar to 
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1 N" 1 N" 1 —5 2" 

^.f h i(a ®a ©a 0a) m ? o h mj a Qa X i X m 



1 , N 2 -N x 
96TT ( -^~> 



and fifteen sums similar to 

1 N " 1 N_1 N_1 

^ i=0 ^ < a 9 a > k=0 **> < a ® a > m=0 V < a G a > X i X k X m 

k^i m^i 

m#c 

1 N(N-l) (N-2) 
" T728 N 3 

Gathering all the terms together, we get, for the 6 moment 



TT m _15_ N(N-l)(N-2) 15 ( N(N-l) v + 1 ( N * 
V a 1728 ^ JJT ; + 960 -^ > * 448 ^3 ; 



V " = 15a 6 + terms in ^ and i- (30) 

w IT 



— 4" — 2" 
Evaluating V V, we set, for example, a = c = d = e, b = f, 

a £ b. Therefore, the sum 
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x N-l 

^ Jo h i, (a b a a © a b) xT 5 " = ^ (-^) 



is the same as before, but many of the other thirty sums are 

zero, leaving only the cases 

i j £ n k m ikjn^m 

i k £ m j n i £ j n m k 

j k t n i m and i m j n k t 

j t m n i k 

j k m n i t 

i j m n k I 

We quickly count 



v 4 „ 2 3 , N(N-1UN-2X , 6 ,N(N-1\ , 1 ,N . 
V a V b " T728 ( ~y } + 960 ( "p^> + 45s ( ^> 



(31) 



V o 4 V K 2 = 3a 6 + terms in ^ and -L 
a b N N 2 



Similarly for V 2 V K 2 V 2 we obtain 
^ a b c 



v 2„ 2„ 2 L_ r N(N-l)(N-2) v 3 r N(N-l) , 1 ,N v 

V a V b v c " T728 ( JJT^ } + 55ff ( ^ > + 55F ( ^T ; 

(32) 

vvv = * 6 + terms in r and 7 
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These results also approach the expected results for 
normal r.v. with an error like 1/N, i.e., 15a , 3a , and a . 

At this point, we comment that the behavior, for large N, 
of any moment is a matter of combinatorial analysis. For example, 
the moment 

V a 2p = 1-3-5--- (2p-l)-a 2p + terms in jjj- , etc. (33) 

can be derived by a counting process. For more complicated 
moments, the analysis is so complicated that we have not under- 
taken it. 

We previously described the N dimensional distribution of 

the W f s in terms of a hypercube which we argued was asymmetically 

N 
distributed among the 2 hyperquadrants in N space. We cannot 

make a simple statement about P (Vq, v p«** v n-1^ ^> ecause 
it is not uniformly distributed, but rather has some density 
which is a continuous function of spacial coordinates. However, 
it is clearly symmetrically distributed in the hyperquadrants, 
by construction. 

Programming Considerations 

The method proposed requires at least the following two or 
three steps: 
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1. Compute N uniform r.v. Xq, X-.,..., X„ -.. Suppose that each 
r.v. requires a time T to generate. Then part 1 will require 
a time NT . 

2. Compute the "Hadamard transform 11 of the array of N numbers 

from part 1. This can be done in-place (i.e., without any extra 

memory) with N log^N additions or subtractions. Assuming that 

an addition and a subtraction take the same amount of time T , 

a' 

part 2 will require a time N(log N)T . If it is desired to 

z a 

adjust the variance of the resulting W , we should add to this 

time an additional time N T , where T is the time for a multiply 

m' m r J 

(or a scale if that is acceptable) . 

3. If the improved approximation to independent gaussian noise 

is desired, generate N additional random variables r. which can 

be used to randomly change or not change the sign of the r.v. 

W. generated in part 2. This should take at most a time 

N(T + T ) assuming the r. can be generated as quickly as a 
u a j 

uniform r.v. and sign changing can be accomplished as quickly 
as addition or subtraction. 

The total time, therefore, is like 



T total =N(2T u + log 2 (2N) T a ) 
but this is the time for generating N r.v. whereas the relevant 
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time would be the time per random variable. This latter time 
is 

T g = 2T U + (log 2 N + 1) T a . 

If many r.v. are needed, as will usually be the case, they 
should be generated N at a time by a subroutine. This sub- 
routine could, of course, supply them one at a time to a user 
program, recomputing a block of N when its supply was exhausted. 
Thus, unless fewer than N r.v. were required, the limitation on 
computing the r.v. N at a time should not be important. 

Conclusions 

Random variables with an almost normal distribution can be 
computed by first obtaining uniform random variables, using 
standard computer algorithms, and then performing a linear trans- 
formation on a set of N = 2 of the uniform random variables to 
obtain N random variables whose joint distribution is not 
gaussian but uniform in a hypercubic region. Nevertheless, many 
of the marginal densities are nearly gaussian, in accord with 
the central limit theorem, and the r.v. are uncorrelated. In 
the limit of large N, all the moments which have been determined 
seem to approach the moments expected for gaussian independent 
r.v. However, some of the joint moments which would be zero for 
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independent gaussian r.v. are finite for the derived set. By 
randomly changing or not changing the sign of the linearly trans- 
formed r.v., a new set of r.v. is obtained for which all the 
moments, including joint moments, show considerable similarity 
to the moments expected from the gaussian independent case, for 
large N. The procedure has the advantage of being faster than 
other decent methods of generating gaussian r.v. on a digital 
computer, with relatively good accuracy. 
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APPENDIX 
NUMBER OF OPERATIONS IN HADAMARD MATRIX MULTIPLICATION 

If we consider the multiplication of an arbitrary vector by 
the matrix H, with elements h. . given by (12), an abvious 
proceedure is to partition the matrix into four quadrants, 



H =■= 






i.j 



i,j 



d. , 






= h. 



iiJ 



b. . 
1>J 




h. . N 
i,j + 2" 


C i,j 


s 


h. , N . 
1 + J, 3 


d U 


s 


h. , N . , N 



N 



* i * f 



- 1 



* j £ I - 1 



It is seen that for 



- ( i^ ( i\ ±l * 1 
a if j - (-1) (-1) 



(-1) 



^-2^-2 



(-1) 



i-l^c-l 



i_i * j _i = for the entire quadrant. Therefore a. . is 



N .. N 



actually an *■ x vr Hadamard matrix. Similarly 



iiJ i-iJ i,j i,J 
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It follows that if the vector to be multiplied by H in a 
column vector [X.] we may form two shorter vectors with components 



t 



X i - X i + X i + f * i * if - 1 



It 



X i -X i -X i+ N * i * | - 1 



N 
and we shall find that the first y components of the product are 



given by 

7 l 

i-0 1 »J x 

N 
and the second -j- components are given by 

N -1 

^ it 

a. . X, 

i-0 1,J x 

each of which are 7 * 7 Hadamard matrix multiplications. This 
reduction cost us N/2 additions and N/2 subtractions. Of course 



a, . may be similarly partitioned and this may be done repeatedly. 

N N 

Each reduction xosts and additional j additions and j subtractions 

so that a total of Nc = N log ? N operations is required in all. 
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